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A new implementation is proposed for including van der Waals interactions in Density Func- 
tional Theory using the Maximally-Localized Wannier functions. With respect to the previous 
DFT/vdW-WF method, the present DFT/vdW-WF2 approach, which is based on the simpler Lon- 
don expression and takes into account the intrafragment overlap of the localized Wannier functions, 
leads to a considerable improvement in the evaluation of the Ce van der Waals coefHcients, as shown 
by the application to a set of selected dimers. Preliminary results on Ar on graphite and Ne on the 
Cu(lll) metal surface suggest that also the C3 coefficients, characterizing molecule-surfaces van der 
Waals interactions are better estimated with the new scheme. 



An accurate description of ubiquitous, long-range van 
der Waals (vdW) interactions is crucial for characteriz- 
ing countless phenomena, belonging to such diverse fields 
as solid state and surface physics, chemistry and biol- 
ogy. For instance, vdW effects are responsible for the 
stabilization of non covalently-bonded crystals and lay- 
ered structures, play a major role in physisorption pro- 
cesses, and are know to affect several biological phenom- 
ena. vdW interactions are due to long range-correlations, 
in particular, the leading term is a consequence of 
correlated, instantaneous dipole fluctuations. Density 
Functional Theory (DFT), thanks to its favorable scal- 
ing properties, represents a popular, efficient and invalu- 
able approach, that is also applicable to extended systems 
where other ab initio schemes turn out to be too compu- 
tationally expensive. However, standard DFT schemes 
only provide a local or semilocal treatment of the elec- 
tronic correlation, so that they are unable to properly 
reproduce genuine vdW effects.— 

The simplest way to include vdW interactions in DFT 
is represented by semiempirical methods,—!^ where, typ- 
ically, an approximately derived Cg/i?"^ term is multi- 
plied by a short-range damping function, with parame- 
ters tailored to the specific system considered. Although 
such an approach is very efficient and often gives a sub- 
stantial improvement with respect to a standard DFT 
method, nonetheless its accuracy is difficult to asses in 
advance and lacks of transferability (for instance, changes 
in atomic polarizabilities by changing the atom environ- 
ment are neglected). Clearly, better reliability, accuracy, 
and transferability can be in principle achieved by adopt- 
ing schemes where vdW corrections are computed by ex- 
ploiting the knowledge of the electronic density distribu- 
tion given by DFT. In recent years several approaches 
have been indeed proposed (for a recent review, see, for 
instance, rcf. 0) In order to circumvent the direct use 
of truly non-local DFT functionals, which are not easy 
to evaluate efficiently, some of these methods introduce 
suitable partitioning schemes into separated interacting 
fragments, either relying on effective atom-atom^^— or 
electronic orbital-orbital^"— pairwise Cq/R^ terms. Al- 
though these techniques are expected to be more reliable 



and transferable than semiempirical approaches, in prac- 
tice, most of them use one or more parameters to be fitted 
using some reference database. 

Here we describe and apply a new implementation 
of the DFT/vdW-WF method^i^iiiii^ where electronic 
charge partitioning is achieved using the Maximally- 
Locahzcd Wannier Functions (MLWFs). The MLWFs 
are obtained from a unitary transformation in the space 
of the occupied Bloch states, by minimizing the total 
spread functionali^^ 
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The localization properties of the MLWFs are of par- 
ticular interest for the implementation of an efficient 
vdW correction scheme: in fact, the MLWFs represent 
a suitable basis set to evaluate orbital-orbital vdW in- 
teraction terms. While in the original DFT/vdW-WF 
method the vdW energy correction for two separate 
fragments was computed using the exchange-correlation 
functional proposed by Andcrsson et al. ^ our novel ver- 
sion (DFT/vdW-WF2 method) is instead based on the 
simpler, well known London's expression>i^ basically, two 
interacting atoms, A and B, are approximated by cou- 
pled harmonic oscillators and the vdW energy is taken 
to be the change of the zero-point energy of the coupled 
oscillations as the atoms approach; if only a single ex- 
citation frequency is associated to each atom, loa, ujb, 
then 
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where Za.b is the total charge of A and B, and Rab 
is the distance between the two atoms (e and m arc the 
electronic charge and mass). Now, adopting a simple 
classical theory of the atomic polarizability, the polariz- 
ability of an electronic shell of charge eZi and mass mZi, 
tied to a heavy undcformable ion can be written as 
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Then, given the direct relation between polarizabihty 
and atomic volume)^ we assume that Ui ^ jSf, where 7 
is a proportionaUty constant, so that the atomic volume 
is expressed in terms of the MLWF spread, Si. Rewrit- 
ing eq. ^ in terms of the quantities defined above, 
one obtains an explicit expression (much simpler than 
the multidimensional integrals involved in the Andersson 
functionalii) for the Cg vdW coefficient: 
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The constant 7 can then be set up by imposing that 
the exact value for the H atom polarizabihty (an =0.866 
a.u.) is obtained (of course, in the H case, one knows the 
exact analytical spread. Si = Sh = a.u.). Note that, 
by expressing the "atomic" volume as a function of Si 
we actually implicitly switch from an atom-atom to an 
orbital-orbital approach. 

In order to achieve a better accuracy, one must prop- 
erly deal with intrafragment MLWF overlap: in fact, the 
DFT/vdW-WF method is strictly valid for nonoverlap- 
ping fragments only; now, while the overlap between the 
MLWFs relative to separated fragments is usually neg- 
ligible for all the fragment separation distances of inter- 
est, the same is not true for the MLWFs belonging to the 
same fragment, which are often characterized by a signifi- 
cant overlap. This overlap affects the effective orbital vol- 
ume, the polarizabihty, and the excitation frequency (see 
eq. ([3])), thus leading to a quantitative effect on the value 
of the Ce coefficient. We take into account the effective 
change in volume due to intrafragment MLWF overlap 
by introducing a suitable reduction factor ^ obtained by 
interpolating between the limiting cases of fully overlap- 
ping and non-overlapping MLWFs. In particular, since 
in the present DFT/vdW-WF2 method the i-th MLWF 
is approximated with a homogeneous charged sphere of 
radius 5"^, then the overlap among neighboring MLWFs 
can be evaluated as the geometrical overlap among neigh- 
boring spheres. To derive the correct volume reweighting 
factor for dealing with overlap effects, we first consider 
the limiting case of two pairs (one for each fragment) 
of completely overlapping MLWFs, which would be, for 
instance, applicable to two interacting He atoms if each 
MLWF just describes the density distribution of a sin- 
gle electron; then we can evaluate a single Cg coefhcicnt 
using eq. (|4]) with Za,b = 2, so that: 
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Alternatively, the same expression can be obtained by 
considering the sum of 4 identical pairwise contributions 



(with Z = 1), by introducing a modification of the effec- 
tive volume in such a way to take the overlap into account 
and make the global interfragment Cg coefficient equiv- 
alent to that in eq. ([5|). This is clearly accomplished by 
replacing Sf in eq. (g]) with £,Sf, where ^ = 1/2. This 
procedure can be easily generalized to multiple overlaps, 
by weighting the overlapping volume with the factor n~^, 
where n is the number of overlapping MLWFs. Finally, 
by extending the approach to partial overlaps, we de- 
fine the free volume of a set of MLWFs belonging to a 
given fragment (in practice three-dimensional integrals 
are evaluated by numerical sums introducing a suitable 
mesh in real space) as: 
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where Wfreei^i) is equal to 1 if |r/ — r^j < S'i for at least 
one of the fragment MLWFs, and is otherwise. 

The corresponding effective volume is instead given by 
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where the new weighting function is defined as Weff{vi) = 
Wfreei't^i) ' 'T-iu(rO~^i with Uwivi) that is equal to the 
number of MLWFs contemporarily satisfying the rela- 
tion |r; — Til < Si. Therefore, the non overlapping por- 
tions of the spheres (in practice the corresponding mesh 
points) will be associated to a weight factor 1, those be- 
longing to two spheres to a 1/2 factor, and, in general, 
those belonging to n spheres to a l/n factor. The average 
ratio between the effective volume and the free volume 
i^eff/^free) IS thcu assigucd to thc factor ^, appearing 
in eq. ([8]). Although in principle the correction factor ^ 
must be evaluated for each MLWF and the calculations 
must be repeated at different fragment-fragment separa- 
tions, our tests show that, in practice, if the fragments 
are rather homogeneous all the ^ factors are very similar, 
and if the spreads of the MLWFs do not change signifi- 
cantly in the range of the interfragment distances of inter- 
est, thc ^'s remain essentially constant; clearly, exploit- 
ing this behavior leads to a significant reduction in the 
computational cost of accounting for thc intrafragment 
overlap. We therefore arrive at the following expression 
for the Cg coefficient: 
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where ^a,b represents the ratio between the effective 
and the free volume associated to the A-ih and B-ih. 
MLWF. The need for a proper treatment of overlap ef- 
fects has been also recently pointed out by Andrinopoulos 
et al.^ who however applied a correction only to very 
closely centred WFCs. 

Finally, the vdW interaction energy is computed as: 
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where f{Rij) is a short-range damping function, which is 
introduced not only to avoid the unphysical divergence 
of the vdW correction at smah fragment separations, but 
also to eliminate double countings of correlation effects 
(in fact standard DFT approaches are able to describe 
short-range correlations); it is defined as: 

The parameter Rs represents the sum of the vdW radii 
= Rf^+Rf^.^ with (by adopting the same criterion 
chosen above for the 7 parameter) 

jyvdW jyvdW /I 1 \ 

" ^ ^ ' 

where R'"jf^ is the literatureii (1.20 A) vdW radius of 
the H atom, and, following Grimme et al.^ a ^ 20 (the 
results are almost independent on the particular value 
of this parameter). Although this damping function in- 
troduces a certain degree of empiricism in the method, 
we stress that a is the only ad-hoc parameter present in 
our approach, while all the others are only determined 
by the basic information given by the MLWFs, namely 
from first principles calculations. 

Calculations were performed using the CPMD codoii^ 
and taking, as the reference DFT GGA functional, both 
the PBEi^ and revPBE^ flavor: PBE is chosen because 
it represents one of the most popular GGA function- 
al for standard DFT simulations of condensed-matter 
systems, while revPBE, which usually gives results close 
to those obtained by a pure Hartree-Fock approach, has 
been used both in our previous DFT/vdW-WF calcula- 
tions and also in other vdW-corrected DFT studiesi^i^ 
As already pointed out elsewhere r^i^iS^ vdW-corrcctcd 
PBE calculations show a general tendency to overbind- 
ing (attributed to an overestimate of the long-range part 
of the exchange contribution), with equilibrium distances 
in reasonable agreement with reference values, while in- 
stead vdW-corrected revPBE typically overestimates the 
equilibrium distances but gives better estimates for the 
binding energies. 

We stress that the computational cost of the 
DFT/vdW-WF2 method, although slightly increased 
with respect to that of the previous DFT/vdW-WF 
scheme, still represents a negligible additional cost if com- 
pared to that of a standard DFT calculation, thus satis- 
fying the basic efficiency requirement. 

In Table H] we report the Cg coefficients computed for 
a set of 18 dimolecular systems, where vdW interactions 
represent the dominant (or at least a significant) contri- 
bution, using our DFT/vdW-WF2 method to be com- 
pared to reference data. As can be seen, in most of the 



systems the Cg coefficient value is reduced and the overall 
performance is much improved with respect to the pre- 
vious DFT/vdW-WF approach)^ in fact the mean rela- 
tive error (MRE) is decreased from 14.1 to 0.3 %, and 
from 16.6 to -0.3 % with revPBE and PBE, respectively, 
while the corresponding reductions in the mean abso- 
lute relative error (MARE) are from 35.4 to 14.6 % with 
revPBE and from 32.8 to 10.8 % with PBE. The effect is 
particularly apparent in rare-gas dimers and dimolecular 
complexes containing benzene, which are systems where 
the correction factor ^ is important due to a significant 
overlap among MLWFs belonging to the same fragment. 
Clearly, one expects that this correction will be also im- 
portant in large molecules and also in extended systems, 
characterized by relatively delocalized electronic charge 
distributions, corresponding to large MLWF spreads. 

Note that, the typical decrease of the Cq coefficient 
values obtained by the DFT/vdW-WF2 method does not 
necessarily lead to a reduction of the vdW energy contri- 
bution; this is clearly due to the effect of the adopted 
damping function, which determines the interplay be- 
tween the Cq/R~^ vdW correction and standard DFT 
energy contributions. As a result, the new DFT/vdW- 
WF2 scheme in general predicts slightly shorter equilib- 
rium distances, in better agreement with reference data 
than DFT/vdW-WF, while instead the binding ener- 
gies exhibit a behavior similar to that obtained by the 
DFT/vdW-WF approach, basically determined by the 
underlying DFT GGA functional (see above comment). 
We also point out that, a much more accurate estimate 
of the Ce coefficients allows for a better description of 
the vdW interactions even for interfragment distances far 
from the equilibrium values, which is of particular rele- 
vance, both for the applications to large systems and for 
Molecular Dynamics simulations. 

In order to test the applicability of the present 
DFT/vdW-WF2 method also to extended systems, which 
of course represent the most interesting application field 
because high-quality chemistry methods are too compu- 
tationally demanding, we considered both the adsorption 
of a single Ar atom on graphite and of a Ne atom on the 
Cu(lll) metal surface, which represent two typical ph- 
ysisorption processes. In the case of Ar on graphite, cal- 
culations have been performed using the same approach 
followed in ref. [l^, while for Ne on Cu(lll) we have used 
the Quantum-ESPRESSO^ ab initio package (MLWFs 
have been generated as a post-processing calculation us- 
ing the WanT package^i): we modeled the substrate us- 
ing a periodically-repeated hexagonal supercell, with a 
(■\/3 X ■\/3)i?30° structure and a surface slab made of 15 
Gu atoms distributed over 5 layers; the Brillouin Zone 
has been sampled using a 6 x 6 x 1 fc-point mesh. 

By fitting the adatom binding energy as a function of 
its distance from the substrate, z, (as it is usually done^^ 
the fit has been performed by optimizing the parame- 
ters of the function: Ae~^^ — ^^/{z — zq)^ ), one can 
easily estimate the C3 coefficients that characterize the 
adatom-surfacc vdW interactions. As can be seen in Ta- 
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ble II, although the agreement with reference C3 data is 
not perfect yet (obtaining accurate C3 coefHcients rep- 
resents a notorious difficult problem, sec, for instance, 
ref. [2^, nonetheless the DFT/vdW-WF2 method gives 
a dramatic improvement with respect to the previous 
DFT/vdW-WF scheme. 

In conclusion, we have described and applied a new 
implementation of our vdW-correction method based 
on the Maximally-Localized Wannier functions: the 
DFT/vdW-WF2 approach is based on the London ex- 



pression and takes into account the MLWF intrafragment 
overlap. The application to selected dimers and also to 
Ar on graphite and Ne on the Cu(lll) metal surface 
show a substantial improvement in the long-range vdW- 
coefBcient (Ce and C3) estimates. Work is in progress to 
achieve a similar level of improvement in equilibrium dis- 
tances and binding energies: this would probably require 
the introduction of more sophisticated, DFT-functional 
dependent, damping functions. 
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DFT/vdW-WF 


DFT/vdW-WF2 


Ref. 


H-H 


7.50(8.0) 


7.17(7.48) 


6.38 


He-He 


0.57(0.62) 


1.48(1.47) 


1.45 


Ne-Ne 


4.35(4.73) 


10.4(8.9) 


6.35 


Ne-Ar 


24.9(16.9) 


26.4(22.9) 


19.5 


Ar-Ar 


92.5(93.2) 


65.8(66.1) 


64.3 


Kr-Kr 


214.0(227.0) 


124.0(124.0) 


131.0 


Xe-Xe 


618.0(621.0) 


261.0(262.0) 


285.9 


N2-N2 


87.4(89.3) 


81.2(80.5) 


73.3 


CO-CO 


85.6(86.7) 


84.8(85.1) 


81.5 


NH3-NH3 


67.1(88.4) 


63.5(77.6) 


89.03 


H20-H20 


35.2(35.6) 


38.9(37.3) 


45.29 


C2H6-C2H6 


308.0(315.0) 


298.0(300.0) 


381.9 


CH4-CH4 


103.0(119.0) 


98.2(111.0) 


129.7 


CsHe-CeHg 


2930.0(2900) 


1710.0(1710) 


1722.7 


CsHs-Ar 


490.0(495.0) 


333.0(334.0) 


330.1 


C6H6-H2O 


323.0(325.0) 


252.0(256.0) 


277.4 


CO2-CO2 


187.0(191.0) 


162.0(158.0) 


158.5 


NH3-CO 


78.5(88.6) 


75.5(80.4) 


90.2 


MRE 


14.1(16.6)% 


0.3(-0.3)% 




MARE 


35.4(32.8)% 


14.6(10.8)% 





TABLE I: Ce coefficients (in meV A''), using the reference 
DFT revPBE functional (PBE in parenthesis) computed with 
the DFT/vdW-WF2 method, compared with those obtained 
by the previous DFT/vdW-WF scheme,— and with reference 



values.- 
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DFT/vdW-WF 


DFT/vdW-WF2 Ref. 


Ar-graphite 


18318 


2057 1210 


Ne-Cu(lll) 


1226 


589 488 



TABLE II: C3 coefficients (in meV A^), using the refer- 
ence DFT revPBE functional computed with the DFT/vdW- 
WF2 method, compared with those obtained by the previous 
DFT/vdW-WF scheme,— and with reference values.— 



